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Abstract 

We present numerical results for a chemical reaction of colloidal particles which are transported 
by a laminar fluid and are focused by periodic obstacles in such a way that the two components 
are well mixed and consequently the chemical reaction is speeded up. The roles of the various 
system parameters (diffusion coefficients, reaction rate, obstacles sizes) are studied. We show that 
focusing speeds up the reaction from the diffusion limited rate ~ i -1 / 2 to very close to the perfect 
mixing rate, ~ t . 
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I. INTRODUCTION 



The motion of colloidal particles on modulated surfaces has attracted a great deal of 

fin 

attention in the past decade [lH6j. The interest in this subject has mainly been directed 
at sorting phenomena, and a considerable portion of the work has been experimental. The 
work has focused on mixtures of particles which are sorted into separate streams of different 
species when the mixture is transported under laminar conditions along a structured or mod- 
ulated medium with periodic obstacles or traps. In this scenario it is possible to control the 
transport of materials such as DNA fragments or functionalized biological colloidal particles. 
The modulated surfaces are specifically designed to present periodic arrays of traps 
or microfabricated obstacles |4[ among other configurations. This technique can be applied 
not only to solid spherical particles but also to other objects such as cells, proteins, DNA, 
and droplets in inmimiscible fluids [7J. Theoretical studies complemented with stochastic 



simulations have received considerable attention |9Hll|. Other sorting methods based on 
inertia and hydrodynamics have also been explored 8|. The sorting phenomenon consists 
of a lateral or orthogonal displacement of the particles with respect to the driving force or 
velocity direction of the fluid mixture. The deviation of the particles in a mixture from 
the direction of flow of the mixture depends on some property or group of properties of 
the particles such as, for instance, size, mass, or charge, causing the particles with different 
values of these properties to flow in different directions. 

The same principle can be used to achieve the converse effect, namely, to focus particles 
coming from different directions and mix them if the modulated structure of obstacles or 

n 

traps is prepared accordingly |5|. Focusing of particles is useful in a number of different 
scenarios such as counting, detecting, and mixing 6[. This property has special relevance 
in the laminar regime, where slow molecular diffusion makes it difficult to concentrate and 
mix particles. 

Our particular focus in this work lies in using this methodology to mix reactants in order 
to speed their reaction. It is our objective to show that using an appropriate modulated or 
structured surface of obstacles one can concentrate two reactants in a very small domain, 
thus favoring their chemical reaction. Our main result is that reactants that arise from 
non-homogeneous distributions can be efficiently mixed and as a result are able to reach 
the classical law of mass action reaction regime characterized by the reactant concentration 
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decay law ~ t^ 1 much sooner than they would in the absence of a mixing mechanism. We 
present results on the efficiency of some obstacle geometries toward this purpose. We discuss 
the roles of the different control parameters and of the densities, the diffusion coefficients of 
the different species, the reaction rates, and the particle and obstacle sizes. 

Our presentation proceeds as follows. First in Sec. [TT] we begin with the description of 
the continuous dynamical scenario and present the associated dynamical equations which 
contain diffusion, advection and reaction terms. In Sec. lIIII we present numerical results from 
the simulation of the equations and thereby analyze the roles of the different parameters. 
Finally we close with some conclusions and perspectives in Sec. [TV] 

II. DYNAMICAL MODEL 

The theoretical scenario is the advective reaction-diffusion model of chemical kinetics 
corresponding to the simple irreversible reaction A + B — > 0. The dynamical equations for 
the two react ants are 

d 

— c a (x, y; t) = -VJ a (x, y; t) - kc a c b , (la) 
d 

— c b (x, y\ t) = -VJ 6 (x, y; t) - kc a c b , (lb) 

where c a and c b are the time dependent local concentrations of the reactants A and B, k is 
the reaction rate constant, and J a , J b are the fluxes of the reactants. The latter are given 
by 

J a (x, y; t) = c a (x, y; t)v(x, y\ t) - D a Vc a - U c a VU, (2a) 
3 b (x, y; t) = c b (x, y; t)v(x, y; t) - D b Vc b - U c b VU. (2b) 

Here D a and D b are the diffusion coeficients, UoU(x/X,y/X) is the modulated potential in- 
teraction due to obstacles, and v(x/A, y/X) is the local velocity responsible for the advective 
flux, which we assume to be a Hele-Shaw flow (that is, a flow between two very close parallel 
plates). We have explicitly extracted the amplitude Uo of the potential so that U(x/\,y/X) 
is the potential of unit amplitude. This potential is modeled by placing a circular tower at 
each obstacle, which changes from unit value at the center of the disk defining the base of 
the obstacle to a zero value outside the range of the interaction. Specifically, it is modeled 
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by the expression 

f(,) = ^(l-tabt^), (3) 

fc=l ^ ' 

where R^, k = 1 . . . N, are the positions of the centers of the N obstacles of base radius 
a, d > a is the radius of the interaction, and 5 is the (small) scale that characterizes a 
substantial change in the value of the potential. To take into account the finite size of the 
advected particles in a model with a continuous concentration field we have introduced an 
interaction potential with a radius d larger than the obstacle radius a. The distance d — a 
then represents the particle radius (see Fig. [TJ. The potential range d corresponds to the 
minimum distance between the centers of the colloidal particle and the obstacle, and hence 
is the sum of their radii. The flow field is obtained by solving the Laplace equation with 
boundary conditions that reflect the presence of the N circular obstacles of radius a, but 
neglecting any effects that the advected particles might have on the flow. Even this is a 
rather arduous task that we have moved to the Appendix. 

Equations ([1]) can be simplified with a change to the new variables r, x', and y', 

t = t r, x = Xx', y = Xy', (4) 

where A and to are characteristic length and time scales. This transforms Eqs. (CQ) to 
d 

— c a (x', y', r) = D a V 2 c a + C/ V(c a W) + V(vc a ) - kc a c b , (5a) 
d 

Q-c b (x', y', r) = D b V 2 c b + U V(c b VU) + V(vc 6 ) - kc a c b . (5b) 
The dimensionless parameters are given by 

and the concentrations and derivative operators are also dimensionless. We fix the parame- 
ters E>i = 0.01 and Uq = 0.01 throughout the paper. From now on for simplicity of notation 
we drop the primes on x' and y' and alert the reader not to confuse x and y as used henceforth 
with the original variables. 

The equations are simulated on a two-dimensional lattice of N x = 1000, N y = 300 cell 
centers and cell dimensions (dimensionless quantities) Ax = Ay = 0.05, which corresponds 
to system dimensions l x = 50 and l y = 15. Our integration time step is At = 5 x 10~ 4 . The 
dynamical evolution is reported every 5 time units up to a final time r = 60. 



FIG. 1: Stream lines of the velocity field (see Appendix), one obstacle (central circle) and the 
motion (arrows) of a particle (small circle) close to the obstacle. The dashed circle indicates the 
area of the obstacle's influence due to the finite size of the particle. The radius of the obstacle a 
and the potential d are indicated. 

Figure [1] shows the role of one obstacle, the flow lines and the finite size of a colloidal 
particle. When the particle following a flux line is close to the obstacle its trajectory changes 
to a different flow line pointing away from the obstacle. The result is a lateral deviation of 
the particle from the initial flow line as it circumnavigates the obstacle. In the presence of 
several obstacles, and depending on their specific spatial distribution, these deviations can 
favor focusing and hence accelerate the reaction. To check this hypothesis we have employed 
two obstacle patterns, one with 293 obstacles in a tilted periodic pattern (PP), as shown 
in Fig. [2J and a second one with the same number of obstacles but distributed randomly 
over the same area (random pattern, RP). We have also performed some simulations with 
no pattern (NP). We have mainly used obstacles of radius a = 0.15 and an interaction 
potential radius d = 0.25, which corresponds to a radius d — a — 0.10 of reacting particles. 
We have also used the values a = 0.10, 0.05 and d = 0.20 in some cases in order to study the 
dependence of the results on these parameters. We have also performed some simulations 
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FIG. 2: Periodic pattern (PP) with tilted obstacles. The black circles are the obstacles and the 
red area corresponds to the range of the interaction potential. The dashed lines show the positions 
where incoming and final flows respectively are numerically measured. The blue and red bars on 
the y axis denote the inlets through which the two reacting especies are introduced into the system 
(see Fig. E]). 

in which the flow is not affected by the obstacles and thus the particles move at a constant 
velocity, i.e., v = constant (a = 0), see below. This is done because in previous calculations 
concerning flows of particle mixtures over surfaces with obstacles we did not take the effect 
of the obstacles on the flow field into account 9Nll| . Here we have the opportunity to assess 
the importance of doing so (see below). 

We are interested in a number of observables along the flow direction, and not in a 
direction perpendicular to the flow, which we integrate over. The observables we focus on 
are the local reaction rate R(x), 

R(x) = k J dyc a (x,y)c h (x } y), (7) 

and the total flux as a function of the position x, 

J(x,t) = J dy[J x>a (x,y,t) + J X)b (x,y;t)] . (8) 

From this quantity the reaction efficiency at a long time (here taken as r = 60) is evaluated 
comparing the fluxes of non reacted particles at two points, x = 0.5 and x = 40, 

J(x = 40) 



FIG. 3: Distribution of the two species, black and white, at the final time r = 60, associated with 
the periodic pattern (PP) of Fig. [2] (top) and with a random pattern (RP) of obstacles (bottom). 
In both cases k = 0.20 . 

In Fig. [2] we present the two-dimensional landscape in which the advection, diffusion, 
and reaction of the two components A and B takes place. One can see the periodic and 
tilted structure of the circular obstacles. Figure [3] (top) shows the distribution of both 
concentrations at time r = 60. It is clear that the obstacles focus the concentrations toward 
the center line, where most of the reaction process takes place. This situation should be 
compared with the case of randomly distributed obstacles, Fig. [3] (bottom). 



III. ANALYSIS OF NUMERICAL RESULTS 



In Fig. H]we present initial and final concentration profiles. The left column shows the 
effect of obstacles on the flow of reactants in the absence of a reaction (k — 0). From top to 
bottom we see the reactant flow patterns when the obstacles are placed in the tilted periodic 
pattern (PP), in the random placement (RP), and without obstacles (NP). We clearly see 
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the strong focusing effect of the PP geometry compared to the other cases. In the right 
column we see the same three cases but now the two species are allowed to react with rate 
coefficient k = 0.2. The reaction clearly proceeds much more rapidly when the reactants 
are focused by the obstacles. The reaction in the presence of random obstacles and of no 
obstacles occurs more slowly and at comparable speeds in the two cases, the former slightly 
more rapidly than the latter. In Fig. [5] the local reaction rate R(x) given in Eq. ([7]) and the 
flux J(x) of Eq. ([9]) in the steady state are plotted for the focusing geometry (PP) along the 
flow direction. We observe repeated positions at which the reaction is more efficient. These 
points correspond to rows where the innermost focusing obstacles are closest together. This 
result opens the possibility of alternative patterns with more active reaction domains. The 
flux is shown for the focusing obstacle geometry for a number of rate coefficients and is seen 
to be a monotonically decaying function of position. For the focusing obstacle geometry 
with k = 0.2, for instance, we calculate a reaction efficiency of ~ 85%, a fairly high value 
for such a relatively small system. In Fig. [6] we explore the dependence of the efficiency 
on the reaction parameter k. We see that the efficiency comes close to a maximum value 
around k = 0.2, with little improvement above that. This almost-independence of the rate 
coefficient is an interesting unanticipated feature. 

It is informative to compare this efficiency with that obtained in the most effective case, 
that is, that of perfect (totally homogeneous) mixing. In this case, assuming equal initial 
concentrations c a = Cb = c. the reaction equation is 

^ = -kc 2 . (10) 

GST 



Its solution is 

c(t) 



(11) 



c(0) 1 + C {0)kr 

We choose the comparison time r = 40. If we further choose c(0) = 1 and k = 0.25, we 
find that c(40)/c(0) = 1/11 which gives an efficiency r] = 0.91. This an upper bound for 
this setup. The focusing configuration thus leads to almost perfect mixing. 

We can use this result to obtain an analytic estimate for the efficiency. If we take the 
velocity of the flowing components ~ 1 we can substitute time for space, r = x/v. We 
expect the flux to have a functional form similar to that of the concentration in Eq. ffTTj) . 



S 




FIG. 4: Left colum: concentration profiles c a (black) and q, (red), at x = 0.5 (dashed lines) and 
x = 40 (solid lines) for r = 60, when the species do not react (k = 0). The upper left panel shows 
the profiles for the periodic pattern (PP), the middle left panel for the random pattern (RP), and 
the lower left panel when there are no obstacles (NP) at all on the surface. Right column: the 
same three cases (PP, RP, and NP)) as in the left column but now incorporating the reaction with 
rate coefficient k = 0.2. 
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X X 

FIG. 5: Left panel: reaction rate parameter R(x) [Eq. (JZJ)] for the configuration PP in the steady 
state (r = 60) for k = 0.2. Right panel: flux [Eq. ([9])] for the same case as left panel, for several 
values of k. From top to bottom: k = 0, 0.025, 0.05, 0.1, 0.2, 0.4. 
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k 

FIG. 6: Reaction efficiency versus reaction rate constant k at the final position xf = 40 (r = 60) 
for the periodic pattern (PP, diamonds), random pattern (RP, circles), and without obstacles (NP, 
triangles). The dashed line corresponds to the reaction efficiency for perfect mixing, cf. Eq. (|lip . 
The full lines are the fittings of Eq. (|12j) . 
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and thus we propose the following expresion for the efficiency at xf- 

v{xF,k) = -^~, (12) 
1 + okxp 

with parameters a and b to be fitted. A nonlinear fit to our numerical data yields a = 1.085 
and b = 1.195 (PP), a = 0.358 and b = 0.788 (RP) and a = 0.105 and b = 0.325 (NP). In 
Fig. Owe have plotted this function for xp = 40. This numerical result implies that most 
of the reacted matter was in the perfect mixing regime characterized by the decay ~ r _1 . 
This is an unexpected result because of the inhomogeneous concentrations in this system. 
It tells us that the reaction is dominated by the small domains near the center of the array 
where there is quasi-perfect mixing, as is seen in Fig. |3j 

To complete our analysis, we will explore the effects on these results of varying the 
diffusion coefficient (assumed to be the same for both species), the focusing pattern geometry, 
and the sizes of obstacles and particles. In the left panel of Fig. [7J we present the effect 
of the diffusion parameter. We indicate the focusing domain for the case of our focusing 
obstacle pattern by dashed lines, that is, the focusing takes place within this domain. We see 
that inside this region the reaction is dominated by focusing, but outside of this region the 
reaction is dominated by diffusion. In the right panel of Fig. [7J we present three fluxes that 
correspond to different geometries. We see that the flux with the focusing pattern decays 
faster than the random pattern which is also more effective than no pattern at all. 

Finally we discuss the effect of the parameters d (radius of potential interaction) and a 
(radius of the obstacles) on the reaction efficiency. From Fig. [8] we can see that the stronger 
effect comes from the change in the parameter d. A comparison of pairs of curves with a 
common value of a shows that the reactant flux is strongly reduced by increasing d from 
0.20 to 0.25. This reduction in the flux implies an increment in the efficiency (see inset 
of Fig. as follows. Returning to Fig. [TJ we see that larger d means that more stream 
lines of the liquid flow are affected by the potential, and as a consequence more particles 
traveling along the stream lines are deviated. Moreover, the deviation of each particle will 
be greater, that is, particles are moved to further stream lines. An increase of d then clearly 
produces a larger focussing effect of the reacting particles, and as a consequence the reaction 
is enhanced. Note that such a change in the parameter d with fixed a corresponds to a change 
in the particle radius d — a (a change in the reactants) while maintaining the same pattern 
of obstacles. 
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FIG. 7: Left panel: plot of the flux for the PP configuration and three values of the diffusion 
coefficient: D = 0.01 (solid line), 0.05 (dashed line) and 0.1 (dotted line). Right panel: fluxes for 
the three configurations considered earlier, from top to bottom: NP, RP and PP (same colors as 
in Fig. [6]). In both panels k = 0.2 and r = 60. 

Next we fix the parameter d and explore the role of the parameter a. In Fig. [8] we see 
that an increase in a reduces the efficiency (see inset). The explanation again lies in the fact 
that the obstacles lead to a deformation of the stream lines (see Fig. [1]), and this effect is 
stronger for larger a. With a stronger deformation the stream lines are displaced further from 
the obstacles, and as a result more particles can escape from the action of the interaction 
potential, whose range d is now held fixed. As a consequence the increase of a at fixed d 
reduces the focussing effect of the pattern and hence the reaction efficiency decreases. 

Perhaps most interesting from an experimental point of view is to analyze the effect 
of varying both a and d while maintaining their difference constant. This corresponds to 
changing the dimensions of the obstacles, easily modified in an experiment, while keeping 
the particle radius d — a constant. This would be the scenario in which one attempts to 
optimize the reaction process for given reactive species. In Fig. [8] we can compare the case 
d = 0.20, a = 0.10 with the case d = 0.25, a = 0.15, which corresponds to changing the 
obstacle size with a common particle radius of d — a = 0.10. Results in the figure show 
that the increase of obstacle size in this example clearly reduces the reactant flux and hence 
improves the reaction efficiency 
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FIG. 8: Comparison of fluxes obtained by using patterns with different radii of obstacles a and of 
potential range d. Solid lines: d = 0.20; dashed lines: d = 0.25. For each d value, different values 
of obstacle radius are represented: from top to bottom a = 0.15 (black), a = 0.10 (red), a = 0.05 
(blue) and a = (green, this case corresponds to the v = const, approximation). Inset: Efficiency 
versus parameter a for d = 0.20 (open symbols) and d = 0.25 (black symbols). 

As a last test, we have compared these results to those of approximation in which the 
deformation of the fluid flow is neglected, i.e., in which the flow velocity is taken to be 
constant and only the scale d of the interaction potential is taken into account. In Fig. [T] 
this situation would correspond to straight, horizontal stream lines. This approximation was 



used in our previous work on particle sorting |9Hll|, and in our present scheme corresponds 
to the limit a = 0. In this situation the stream lines are not deformed and hence a larger 
number of particles are deviated by the potential. In Fig. [8] we see the expected result, 
namely, that this approximation overestimates the reaction efficiency when compared to a 
finite a value (for a given value of d). For the smallest values of a (obstacles smaller than 
the particles) this approximation does not qualitatively change the focussing scenario, but 
as a increases the importance of including the effects of advection clearly increases. 
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IV. CONCLUSIONS 



Obstacles placed in carefully selected geometrical patterns have been used in the past to 
effectively separate colloidal mixtures that are caused to flow over a surface containing such 
obstacles. In this paper we have explored the converse, namely, the possibility of speeding 
up the mixing of components that undergo advective diffusion over a surface containing 
carefully situated obstacles. We have illustrated the effects of this mixing by considering 
the reaction of two species and comparing the reaction rates when the species are allowed 
to mix by ordinary advective diffusion in the absence and presence of these obstacles. We 
have shown that a periodic pattern of tilted obstacles, as opposed to a random placement, 
is able to effectively focus the streams of reactive species. We have furthermore shown that 
this focusing mechanism leads rather rapidly to reaction rates comparable to those obtained 
with perfect mixing. We have studied the dependence of the reaction efficiency on the 
different parameters of the problem. We interpret the focusing mechanism as a consequence 
of the finite size of the reacting particles. These results could be useful for the design and 
interpretation of new experiments. 

This work is supported by Ministerio de Economia y Competitividad (Spain) through 
project FIS2012-37655, by the Generalitat de Catalunya through project 2009SGR-878, and 
by the US National Science Foundation through Grant No. PHY-0855471. 



Appendix: Numerical solution of the Hele-Shaw problem with circular obstacles 

We will consider a fluid contained between two parallel plates, separated by a gap d (Hele- 
Shaw cell). In the limit of small d the velocity field u can be considered as two-dimensional, 
and is given by 

d 2 

u = -— Vp, (A.l) 

which is identical to the Darcy's law for the flow in a porous medium. \i is the viscosity of 
the fluid and p is the pressure, which satisfies V 2 p = 0. We place N circular obstacles or 
disks of radius a centered at positions R fe , k — 1 . . . N. The velocity of the fluid far from the 
obstacles is Uqo. At the rigid boundaries the velocity satisfies the condition of zero normal 
component, but not the no-slip boundary condition. The stream lines passing the obstacles 
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are identical to those of a two dimensional inviscid fluid with the same geometry 12]. 
The velocity field can be written in terms of the velocity potential <j) as 



u = V0, (A.2) 

which obeys the Laplace equation 

V 2 = 0. (A.3) 
The general solution of the 2D Laplace equation ( ]A.3j) can be written in polar coordinates 



as 

oo 



/ A — A\ 

A=l ^ a a / 

where axi, b\i are constants that depend on the boundary conditions. It is easy to establish 
that for real these constants should satisfy 

axibxi = (axib\ 2 )* (A. 5) 

CJA2&A2 = (aA2^Al)*- (A. 6) 

We will use the multipole expansion ( 1A.4I) for the velocity field in the neighbhood of each 
disk, with the polar coordinates centered at the disk. Hence we will have N such expansions. 
The advantage of these expansions is that the boundary conditions at the disks are easy to 
formulate. In particular, the normal velocity vanishes, that is, 

= 0, (A.7) 



r=a 



dr 

and then, taking into account the conditions Eqs. (1A.5j) and (1A.6j) we get 



a-xi = a\2 = 1 (A. 8) 

&A1 = b* X2 EE & A . (A.9) 

We still have to find one set of 6a constants for each disk k (k = 1 . . .N), which we will 
denote as b\(k). 

We next consider the ensemble of A^ disks. All multipole expansions should simultane- 
ously satisfy boundary conditions on all the discs and at infinity. It is convenient to employ 
a complex variable for the position, i.e., we define the complex position f as 

r — x + yi. (A. 10) 
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Henceforth a tilde on any vector will denote a similar definition as a complex variable. By 
employing this notation the general solution Eq. flA.4j) can be expressed more compactly as 
a series of powers of f. 

Note that each expansion ()A.4j) is valid near the corresponding disk, but not far from 
it. In particular the expansions diverge at infinity. The way to manage an expression 
for the potential valid at arbitrary distances is to use only the decreasing powers of the 
multipole expansions, adding the terms corresponding to all the disks. Taking into account 
the boundary condition at infinity we can then write 

N 

0(r) = r- Uoo + ^0 fc (r), (A.ll) 

k=l 

where 0fc(r) is the distortion of the potential produced by the disc k, which will have the 
form 

oo r , v \ 

" " (A.12) 



»*(r) = E 

A=l 



c\(k) (^-] + c.c. 



Here f*. = r — is the position relative to the center of disk k, and the c\(k) are constants. 
With this expression the boundary condition at infinity has already been taken into account. 
To find the unknowns c\(k) we should to make this expansion compatible with those of each 
disk, Eq. (1A.4I) . which incorporating boundary contitions can be written as 

= E ( b ^) (jr) " + KU) (7) " + o.cA , j = 1 . . . N. (A.13) 



To relate the expansion Eq. (lA.llj) with Eq. (IA.12j) to those of Eq. flA.13j) . it is convenient 



to use the Taylor expansion of the negative powers of f& (position relative to disk k) in terms 
of the positive powers of fj (position relative to any other disk j 7^ k). The idea is that 
the positive powers in the local expansion Eq. ( ]A.13j) of any disk should correspond to the 
terms coming from the rest of the disks in Eq. (1A.11I) . That is, we write 



y N / \ A+A' 




where R^j = Rj — Rk so that = R^ + fj, and the constants ^a i a can be obtained from 
the recurrence relation 



<?A'A = qy x-i 



1 - A - A' 
A 

qx>o = I- (A.15) 
16 



We now substitute the expansion Eq. (1A.14|) into Eq. (1A.12j) . to obtain 

A 



1 „ f 
( r ) = 2 aU °°a 



+ E 



A=l 



EE 

fc^j A' = l 



+ C.C. 



A+A' 




) 


(?) 







A=0 

By equating terms of the same power of f in Eqs. HA.13[) and (1A. 16|) we find 
bxU) = c a(j), A = 1...JV, 



j = l...N. (A.16) 



A' 



k^j A'=l 



(I 



a „ 
—it 



Rkj 



A'+l 



fc^i A'=l 

oo 



a 

A'+A 



A 



a 



Rki 



A > 1. 



(A.17) 
(A.18) 

(A.19) 

(A.20) 



k^j X'=l \- l kj 

These equation can be solved to obtain the constants c\(J) with which the potential 
and the velocity field can be calculated through Eq. OA. lip . In practice one can use a few 
multipole terms, A = 1 . . . X ma x, for a moderate value of \ ma x since the series converge rapidly 
(more terms are necessary as the distances between disks are decreased). The procedure 
involves finding c\(j) only once, and permits us to find fluid velocity at any position by only 
summing a few contributions from each disk. 

A dipolar approximation {\ max = 1) provides reasonably good results if the disks are 
not very close to each other, and its solution can be used as an initial step of an iterative 
solution of the complete problem. This solution reads: 

2 



dO^ffioo-^oE hp: i a = i....v. 

k^j \ k , 



(A.21) 



Then the velocity field in this approximation is 



N 



I —35- + C.C. 



A' 



1 3 



u y (r) ~ u yoo - a ^ ( + c.c. J 

7 = 1 ^ r 3 ' 



(A.22) 
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Finally, the system of equations (1A.18j) -( fA.20j) takes the following form: 

/ \ A ' +1 



a 



k+j A'=l 



R kj. 

A'+A 



CA 



U) = EE c ^) 



q\'X 



k^j A'=l 

With these, the velocity field is given by 

N oo 



K kj, 



A > 1. 



a 



A 



«,(r) = Ux oo "EE p A(j) ( r ~_^.)A + i 



+ c.c 



AT oo 

%(r) = u yoo - E ( Aca( ^7^ 
i=i A'=i 



r/ 



A 



i + c.c 



(A.23) 
(A.24) 



(A.25) 



To solve Eqs. ( 1A.23|) and ( 1A.24j) an iterative method can be used. These equations can 
be formulated ClS db matrix relation, 



C = AC* + B, 



(A.26) 



where C is a vector containing all the c\(k) coefficients. Then we can apply the following 
iteration, which can be shown to be convergent: 



Q = ACU + B, 



(A.27) 



with the initial term Co given by the dipolar approximation Eq. ( 1A.21I) 
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